In [1]:
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')
In [5]:
x = np.load('./data.npy')
Fs = 512.
In [7]:
plt.figure(figsize=(15,5))
plt.plot(np.arange(0,10+1/Fs,1/Fs),x,'k')
plt.xlabel('time (s)')
Out[7]:
In [46]:
# Process signal
from scipy.signal import firwin, filtfilt
def mybandpass(x,cf,Fs, w = 3):
numtaps = w * Fs / np.float(cf[0])
if numtaps % 2 == 0:
numtaps = numtaps + 1
taps = firwin(numtaps, np.array(cf) / np.float(Fs) * 2, pass_zero=False)
return filtfilt(taps,[1],x)
cf = (4,200)
ecogBhi = mybandpass(ecogB,cf,Fs)
ecogDhi = mybandpass(ecogD,cf,Fs)
In [52]:
from scipy.signal import butter
def _butterstimrmv(x, cf, bw, rate=1000, order=3):
'''
Notch Filter the time series x with a butterworth with center frequency cf
and bandwidth bw
'''
nyq_rate = rate / 2.0
f_range = [cf - bw / 2.0, cf + bw / 2.0]
Wn = (f_range[0] / nyq_rate, f_range[1] / nyq_rate)
b, a = butter(order, Wn, 'bandstop')
return filtfilt(b, a, x)
In [213]:
stimfreq = [65,128,195]
stimbw = [2,30,1]
Nfilters = len(stimfreq)
for nf in range(Nfilters):
if nf == 0:
ecogB2 = _butterstimrmv(ecogBhi,stimfreq[nf],stimbw[nf],rate=Fs)
ecogD2 = _butterstimrmv(ecogDhi,stimfreq[nf],stimbw[nf],rate=Fs)
else:
ecogB2 = _butterstimrmv(ecogB2,stimfreq[nf],stimbw[nf],rate=Fs)
ecogD2 = _butterstimrmv(ecogD2,stimfreq[nf],stimbw[nf],rate=Fs)
In [214]:
from voytoys.scott.spec_tools import fftmed
Hzmed = 0
xlim = (0,205)
f, psdB = fftmed(ecogB, Fs=Fs, Hzmed = Hzmed)
f, psdD = fftmed(ecogD, Fs=Fs, Hzmed = Hzmed)
f, psdBhi = fftmed(ecogBhi, Fs=Fs, Hzmed = Hzmed)
f, psdDhi = fftmed(ecogDhi, Fs=Fs, Hzmed = Hzmed)
f, psdB2 = fftmed(ecogB2, Fs=Fs, Hzmed = Hzmed)
f, psdD2 = fftmed(ecogD2, Fs=Fs, Hzmed = Hzmed)
In [215]:
plt.figure(figsize=(15,10))
plt.subplot(3,1,1)
plt.plot(f,np.log10(psdB),'k',label='before')
plt.plot(f,np.log10(psdD),'r',label='DBS')
plt.legend(loc='best')
plt.ylabel('Raw')
plt.xlim(xlim)
plt.subplot(3,1,2)
plt.plot(f,np.log10(psdBhi),'k',label='before')
plt.plot(f,np.log10(psdDhi),'r',label='DBS')
plt.ylabel('Bandpass 4-200Hz')
plt.xlim(xlim)
plt.subplot(3,1,3)
plt.plot(f,np.log10(psdB2),'k',label='before')
plt.plot(f,np.log10(psdD2),'r',label='DBS')
plt.ylabel('High frequency\n artifacts removed')
plt.xlim(xlim)
Out[215]:
In [229]:
Hzmed = 1
f, psdB2 = fftmed(ecogB2, Fs=Fs, Hzmed = Hzmed)
f, psdD2 = fftmed(ecogD2, Fs=Fs, Hzmed = Hzmed)
plt.figure(figsize=(15,6))
plt.plot(f,np.log10(psdB2),'k',label='before')
plt.plot(f,np.log10(psdD2),'r',label='DBS')
plt.ylabel('High frequency\n artifacts removed')
plt.xlim(xlim)
Out[229]:
In [216]:
t = np.arange(0,10+1/Fs,1/Fs)
trange = (0,1)
tidx = np.logical_and(t>=trange[0],t<trange[1])
plt.figure(figsize=(15,5))
plt.subplot(2,1,1)
plt.plot(t[tidx],ecogB[tidx]-np.mean(ecogB[tidx]),'k',label='raw')
plt.plot(t[tidx],ecogBhi[tidx],'r',label='4-200Hz')
plt.plot(t[tidx],ecogB2[tidx],'b',label='artifacts\nremoved')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(t[tidx],ecogD[tidx]-np.mean(ecogD[tidx]),'k',label='raw')
plt.plot(t[tidx],ecogDhi[tidx],'r',label='4-200Hz')
plt.plot(t[tidx],ecogD2[tidx],'b',label='artifacts\nremoved')
Out[216]:
In [217]:
from pacpy.pac import comodulogram
from matplotlib import cm
# Parameters
fp = (6,40)
fa = (20,200)
dp = 4
da = 8
pac_method = 'mi_tort'
f_phases = np.arange(fp[0], fp[1], dp)
f_amps = np.arange(fa[0], fa[1], da)
In [218]:
# Calculate comods
comodB = comodulogram(ecogB, ecogB, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodD = comodulogram(ecogD, ecogD, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodBhi = comodulogram(ecogBhi, ecogBhi, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodDhi = comodulogram(ecogDhi, ecogDhi, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodB2 = comodulogram(ecogB2, ecogB2, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodD2 = comodulogram(ecogD2, ecogD2, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
In [219]:
climraw = (0,.01)
climhi = (0,.01)
clim2 = (0,.01)
plt.figure(figsize=(15,15))
plt.subplot(3, 3, 1)
plt.pcolor(f_phases, f_amps+da, comodB.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.title('before DBS', size=20)
plt.yticks(np.arange(50,250,50),size=20)
plt.xticks(np.arange(10,40,10),size=20)
plt.xlabel('Phase frequency (Hz)', size=20)
plt.ylabel('Amplitude frequency (Hz)', size=20)
plt.subplot(3, 3, 2)
plt.pcolor(f_phases, f_amps+da, comodD.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.title('DBS', size=20)
plt.subplot(3, 3, 3)
plt.pcolor(f_phases, f_amps+da, comodD.T - comodB.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.title('change', size=20)
plt.subplot(3, 3, 4)
plt.pcolor(f_phases, f_amps+da, comodBhi.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.subplot(3, 3, 5)
plt.pcolor(f_phases, f_amps+da, comodDhi.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.subplot(3, 3, 6)
plt.pcolor(f_phases, f_amps+da, comodDhi.T - comodBhi.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.subplot(3, 3, 7)
plt.pcolor(f_phases, f_amps+da, comodB2.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.subplot(3, 3, 8)
plt.pcolor(f_phases, f_amps+da, comodD2.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.subplot(3, 3, 9)
plt.pcolor(f_phases, f_amps+da, comodD2.T - comodB2.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.tight_layout()
figname = 'Charlescomod'
plt.savefig('C:/Users/Scott/Google Drive/Voytek/Voytek Lab Shared resources/manuscripts/PD/feb2016fig/'+figname + '.png')
In [220]:
from pacpy.pac import pa_series
f_lo = (13,30)
f_hi = (50,200)
phaB2, ampB2 = pa_series(ecogB2, ecogB2, f_lo, f_hi, fs=Fs)
phaD2, ampD2 = pa_series(ecogD2, ecogD2, f_lo, f_hi, fs=Fs)
def _edgeadd_paseries(amp, fosc, Fs, w = 3):
"""
Undo the removal of edge artifacts done by pacpy in order to align
the extrema with their amplitudes
"""
Ntaps = np.floor(w * Fs / fosc[0])
amp2 = np.zeros(len(amp)+2*Ntaps)
amp2[Ntaps:-Ntaps] = amp
return amp2
phaB2 = _edgeadd_paseries(phaB2,f_lo,Fs)
ampB2 = _edgeadd_paseries(ampB2,f_lo,Fs)
phaD2 = _edgeadd_paseries(phaD2,f_lo,Fs)
ampD2 = _edgeadd_paseries(ampD2,f_lo,Fs)
In [221]:
trange = (0,5)
tidx = np.logical_and(t>=trange[0],t<trange[1])
plt.figure(figsize=(15,6))
plt.subplot(3,1,1)
plt.plot(t[tidx],ecogB2[tidx])
plt.ylabel('Raw voltage')
plt.subplot(3,1,2)
plt.plot(t[tidx],ampB2[tidx])
plt.ylabel('HFA amp')
plt.subplot(3,1,3)
plt.plot(t[tidx],phaB2[tidx])
plt.ylabel('Beta phase')
plt.figure(figsize=(15,6))
plt.subplot(3,1,1)
plt.plot(t[tidx],ecogD2[tidx])
plt.ylabel('Raw voltage')
plt.subplot(3,1,2)
plt.plot(t[tidx],ampD2[tidx])
plt.ylabel('HFA amp')
plt.subplot(3,1,3)
plt.plot(t[tidx],phaD2[tidx])
plt.ylabel('Beta phase')
Out[221]:
In [245]:
# Simulate a low and high frequency oscillation.
def norm01(x):
return (x - np.min(x))/(np.max(x)-np.min(x))
from voytoys.scott.tools_sim import simphase
T = 2
Fs = 1000.
t = np.arange(0,T,1/Fs)
betaphase, beta = simphase(T, f_lo, dt=1/Fs, returnwave=True)
_, gamma = simphase(T, f_hi, dt=1/Fs, returnwave=True)
pacsig = norm01(beta) + gamma*norm01(np.cos(betaphase))*.2
plt.figure(figsize=(15,5))
plt.plot(t,pacsig)
plt.xlim((0,1))
fpac,psdpac = fftmed(pacsig, Fs=Fs, Hzmed = 0)
plt.figure()
plt.plot(fpac,np.log10(psdpac))
plt.xlim((0,250))
Out[245]:
In [246]:
comodpac = comodulogram(pacsig,pacsig, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
In [248]:
plt.figure(figsize=(7,7))
plt.pcolor(f_phases, f_amps+da, comodpac.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
plt.yticks(np.arange(50,250,50),size=20)
plt.xticks(np.arange(10,40,10),size=20)
plt.xlabel('Phase frequency (Hz)', size=20)
plt.ylabel('Amplitude frequency (Hz)', size=20)
Out[248]:
In [104]:
# Calculate PAC
from pacpy.pac import ozkurt
pacB = ozkurt(ecogB2, ecogB2, f_lo, f_hi, fs=Fs)
pacD = ozkurt(ecogD2, ecogD2, f_lo, f_hi, fs=Fs)
In [119]:
print pacB, pacD
In [111]:
# Identify peaks and troughs
from voytoys.shape.shape import findpt, Esharp, ESRsharp, RDsteep
boundaryS = 100
pkB, trB = findpt(ecogB2, f_lo, f_lopass = None, Fs = Fs, boundary = boundaryS)
pkD, trD = findpt(ecogD2, f_lo, f_lopass = None, Fs = Fs, boundary = boundaryS)
In [114]:
# Calculate sharpness of extrema
ampPC = 0
spPB = Esharp(ecogB2, pkB, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)
spTB = Esharp(ecogB2, trB, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)
spPD = Esharp(ecogD2, pkD, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)
spTD = Esharp(ecogD2, trD, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)
In [109]:
# Calculate extrema sharpness ratio on multiple time scales
esrmethod = 'aggregate'
width = 5
esrB = np.log10(ESRsharp(ecogB2,pkB,trB,width,esrmethod=esrmethod))
esrD = np.log10(ESRsharp(ecogD2,pkD,trD,width,esrmethod=esrmethod))
In [110]:
# Calculate rise & decay steepness
def RDsteepratio(Rsteep,Dsteep,logrds = True):
if logrds:
return np.log10(np.max((np.mean(Rsteep)/np.mean(Dsteep),np.mean(Dsteep)/np.mean(Rsteep))))
else:
return np.max((np.mean(Rsteep)/np.mean(Dsteep),np.mean(Dsteep)/np.mean(Rsteep)))
RsteepB, DsteepB = RDsteep(ecogB2, pkB, trB)
RDsteeprB = RDsteepratio(RsteepB,DsteepB)
RsteepD, DsteepD = RDsteep(ecogD2, pkD, trD)
RDsteeprD = RDsteepratio(RsteepD,DsteepD)
In [118]:
plt.figure(figsize=(5,8))
plt.subplot(2,1,1)
Nbins=np.arange(-2,2,.1)
plt.hist(np.log10(spPB),Nbins,color='k',label='peaks',alpha=0.5)
plt.hist(np.log10(spTB),Nbins,color='r',label='troughs',alpha=0.5)
plt.ylabel('Count',size=20)
plt.legend(loc='best')
ax1 = plt.gca()
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax1.get_yticklabels(), visible=False)
plt.subplot(2,1,2)
plt.hist(np.log10(spPD),Nbins,color='k',label='peaks',alpha=0.5)
plt.hist(np.log10(spTD),Nbins,color='r',label='troughs',alpha=0.5)
plt.xlabel('Log Extrema sharpness',size=20)
plt.xticks(np.arange(-1,2,1),size=20)
plt.ylabel('Count',size=20)
plt.legend(loc='best')
ax1 = plt.gca()
plt.setp(ax1.get_yticklabels(), visible=False)
print sp.stats.ttest_ind(spPB,spTB)
print sp.stats.ttest_ind(spPD,spTD)
In [ ]: